// eee.Sheffield.PZ.Imaging
//
// Copyright ?Ping Zou, 2007
// sg71.cherub@gmail.com

namespace eee.Sheffield.PZ.Imaging
{
    using System;
    using System.IO;
    using System.Drawing;
    using System.Drawing.Imaging;
    using eee.Sheffield.PZ.Math;
    using System.Collections.Generic;
    using eee.Sheffield.PZ.Imaging.ImageFilter;

    public class NetworkJoint
    {
        public List<TrackResult> _parents = new List<TrackResult>();
        public List<TrackResult> _children = new List<TrackResult>();
    }

    public class RootToLeafFlowRoute
    {
        public List<TrackResult> _nodes = new List<TrackResult>();
        public void Add(TrackResult tr)
        {
            _nodes.Add(tr);
        }
    }

    /// <summary>
    /// post process tracking results
    /// objects:
    /// 1. render
    /// 2. clear duplicates/too shorts ones/wrong detections
    /// 3. connect vessel segments
    /// </summary>
    public class PostProcessingTrackingResults
    {
        #region Field
        private List<TrackResult> _srcTrackResultList;
        private List<TrackResult> _dstTrackResultList;

        private Bitmap _srcImage;
        private Bitmap _dstImage;

        private string _pathName = "";
        private int _lengthThreshold = 3;

        private List<NetworkJoint> _networkJointList = new List<NetworkJoint>();
        private List<RootToLeafFlowRoute> _rootToLeafFlowRouteList = new List<RootToLeafFlowRoute>();
        #endregion

        #region Properties
        public List<TrackResult> SrcTrackResultList
        {
            get { return _srcTrackResultList; }
            set { _srcTrackResultList = value; }
        }
        public List<TrackResult> DstTrackResultList
        {
            get { return _dstTrackResultList; }
            set { _dstTrackResultList = value; }
        }
        public Bitmap SrcImage { get { return _srcImage; } set { _srcImage = value; } }
        public Bitmap DstImage { get { return _dstImage; } set { _dstImage = value; } }
        public string PathName { get { return _pathName; } set { _pathName = value; } }
        public int LengthThreshold { get { return _lengthThreshold; } set { _lengthThreshold = value; } }
        #endregion

        #region Constructor
        /// <summary>
        /// constructor
        /// </summary>
        /// <param name="trackResultList"></param>
        /// <param name="srcImage"></param>
        public PostProcessingTrackingResults(List<TrackResult> trackResultList, Bitmap srcImage)
        {
            _srcTrackResultList = new List<TrackResult>(trackResultList);
            _dstTrackResultList = new List<TrackResult>(trackResultList);
            _srcImage = new Bitmap(srcImage);
            eee.Sheffield.PZ.Imaging.Image.StandardizeImageTo24bppRGB(ref _srcImage);
            _dstImage = new Bitmap(_srcImage);
            if (_dstImage.PixelFormat == PixelFormat.Format24bppRgb)
                eee.Sheffield.PZ.Imaging.Image.StandardizeImageTo24bppRGB(ref _dstImage);
        }

        public PostProcessingTrackingResults()
        {

        }
        #endregion

        #region post process methods

        public bool CheckRoutes()
        {
            return CheckRoutes(_rootToLeafFlowRouteList, _dstTrackResultList);
        } // CheckRoutes()

        /// <summary>
        /// check if routes contain all track results
        /// </summary>
        /// <param name="rootToLeafFlowRouteList"></param>
        /// <param name="trackResultList"></param>
        /// <returns></returns>
        private bool CheckRoutes(List<RootToLeafFlowRoute> rootToLeafFlowRouteList, List<TrackResult> trackResultList)
        {
            bool searched = false;
            int currentIndex;
            foreach (TrackResult tr in trackResultList)
            {
                searched = false;
                currentIndex = tr._index;
                foreach (RootToLeafFlowRoute route in rootToLeafFlowRouteList)
                {
                    foreach (TrackResult tr1 in route._nodes)
                        if (tr1._index == currentIndex)
                        {
                            searched = true;
                            break;
                        }
                    if (searched)
                        break;
                }
                if (!searched)
                    break;
            }

            return searched;
        } // CheckRoutes()

        public void FindRootToLeafRoute()
        {
            _rootToLeafFlowRouteList = FindRootToLeafRoute(_dstTrackResultList);
        } // FindRootToLeafRoute()

        /// <summary>
        /// find roots to leave routes
        /// </summary>
        /// <param name="trackResultList"></param>
        /// <returns></returns>
        private List<RootToLeafFlowRoute> FindRootToLeafRoute(List<TrackResult> trackResultList)
        {            
            List<RootToLeafFlowRoute> rootToLeafFlowRoute = new List<RootToLeafFlowRoute>();

            int totalRoute = trackResultList.Count;
            
            int currentRouteCount;
            List<RootToLeafFlowRoute> routeToSearch = new List<RootToLeafFlowRoute>();
            RootToLeafFlowRoute routeToAdd;
            RootToLeafFlowRoute currentRoute;
            TrackResult currentNode;
            List<TrackResult> childrenNodeList = new List<TrackResult>();
            // start from roots
            foreach (TrackResult root in trackResultList)
            {
                if (root._isRoot)
                {
                    routeToAdd = new RootToLeafFlowRoute();
                    routeToAdd.Add(root);
                    routeToSearch.Add(routeToAdd);
                    
                    while (routeToSearch.Count > 0)
                    {
                        currentRouteCount = routeToSearch.Count;
                        // search downstream each route
                        for (int routeIndex = 0; routeIndex < currentRouteCount; routeIndex++)
                        {
                            currentRoute = routeToSearch[routeIndex];
                            currentNode = currentRoute._nodes[currentRoute._nodes.Count - 1];
                            if (currentNode._isLeaf)
                            {
                                // terminate the current route
                                rootToLeafFlowRoute.Add(currentRoute);
                                routeToSearch.Remove(currentRoute);
                                totalRoute--;
                            }
                            else
                            {
                                // get current node's children
                                foreach (TrackResult tr in trackResultList)
                                    if (tr._rootIndexList.Contains(currentNode._index))
                                        childrenNodeList.Add(tr);                                
                                // add routes
                                for (int i = 1; i < childrenNodeList.Count; i++)
                                {
                                    routeToAdd = new RootToLeafFlowRoute();
                                    routeToAdd._nodes.AddRange(currentRoute._nodes);
                                    routeToAdd.Add(childrenNodeList[i]);
                                    routeToSearch.Add(routeToAdd);
                                }
                                currentRoute.Add(childrenNodeList[0]);
                            }

                            if (totalRoute == 0)
                                break;
                        } // for (int routeIndex = 0; routeIndex < currentRouteCount; routeIndex++)
                        if (totalRoute == 0)
                            break;
                    } // while (routeToSearch.Count > 0)
                } // if (root._isRoot)
                if (totalRoute == 0)
                    break;
            } // foreach (TrackResult root in trackResultList)

            return rootToLeafFlowRoute;
        } // FindRootToLeafRoute()
        
        public bool CheckJoint()
        {
            return CheckJoint(_networkJointList, _dstTrackResultList);
        } // CheckJoint()

        /// <summary>
        /// check network joint list contains all track results
        /// </summary>
        /// <param name="networkJointList"></param>
        /// <param name="trackResultList"></param>
        /// <returns></returns>
        private bool CheckJoint(List<NetworkJoint> networkJointList, List<TrackResult> trackResultList)
        {
            bool searched = false;
            int currentIndex;
            foreach (TrackResult tr in trackResultList)
            {
                searched = false;
                currentIndex = tr._index;
                foreach (NetworkJoint joint in networkJointList)
                {
                    foreach (TrackResult parent in joint._parents)
                        if (parent._index == currentIndex)
                        {
                            searched = true;
                            break;
                        }
                    foreach (TrackResult child in joint._children)
                        if (child._index == currentIndex)
                        {
                            searched = true;
                            break;
                        }
                    if (searched)
                        break;
                }

                if (!searched)
                    break;
            }

            return searched;
        } // CheckJoint()

        public bool CheckJointFlow()
        {
            return CheckJointFlow(_networkJointList);
        } // CheckJointFlow()

        /// <summary>
        /// check joint flow balance, in flow = out flow
        /// </summary>
        /// <param name="networkJointList"></param>
        /// <returns></returns>
        private bool CheckJointFlow(List<NetworkJoint> networkJointList)
        {
            double jointInFlow;
            double jointOutFlow;
            bool isBalance = true;
            foreach (NetworkJoint joint in networkJointList)
            {
                jointInFlow = 0.0;
                jointOutFlow = 0.0;
                foreach (TrackResult tr in joint._parents)
                    jointInFlow += tr._flow;
                foreach (TrackResult tr in joint._children)
                    jointOutFlow += tr._flow;
                if (jointInFlow != jointOutFlow)
                {
                    isBalance = false;
                    break;
                }
            }

            return isBalance;
        } // CheckJointFlow()
        
        public void InitializeFlow()
        {
            _networkJointList = InitializeFlow(_dstTrackResultList);
        } // InitializeFlow()

        /// <summary>
        /// initialize flow to each pipe
        /// </summary>
        /// <param name="trackResultList"></param>
        /// <returns></returns>
        private List<NetworkJoint> InitializeFlow(List<TrackResult> trackResultList)
        {
            List<NetworkJoint> networkJointList = new List<NetworkJoint>();
            TrackResult current;
            int currentIndex;
            double jointInFlowSum;
            double jointOutFlowSum;
            List<TrackResult> children = new List<TrackResult>();
            List<TrackResult> parents = new List<TrackResult>();

            List<TrackResult> downstreamSearch = new List<TrackResult>();
            List<TrackResult> upstreamSearch = new List<TrackResult>();
            int unknownCount;
            bool isIn;

            // start from #1
            current = trackResultList[0];
            if (current._index != 1)
            {
                foreach (TrackResult tr in trackResultList)
                {
                    if (tr._index == 1)
                    {
                        current = tr;
                        break;
                    }
                }
            }
            downstreamSearch.Add(current);

            while (true)
            {
                // downstream search
                if (downstreamSearch.Count > 0)
                {
                    current = downstreamSearch[0];
                    currentIndex = current._index;
                    // get current's children
                    children.Clear();
                    foreach (TrackResult tr in trackResultList)
                    {
                        foreach (int rootIndex in tr._rootIndexList)
                        {
                            if (rootIndex == currentIndex && !children.Contains(tr))
                            {
                                children.Add(tr);
                            }
                        }
                    }
                    
                    // get "childern"'s parents
                    parents.Clear();
                    parents.Add(current);
                    foreach (TrackResult child in children)
                    {
                        foreach (int rootIndex in child._rootIndexList)
                        {
                            if (rootIndex != currentIndex)
                            {
                                foreach (TrackResult tr in trackResultList)
                                {
                                    if (tr._index == rootIndex && !parents.Contains(tr))
                                    {
                                        parents.Add(tr);
                                        break;
                                    }
                                }
                            }
                        }
                    }

                    // assign flow, downstream, => start from parent
                    jointInFlowSum = 0.0;
                    jointOutFlowSum = 0.0;
                    unknownCount = 0;
                    foreach (TrackResult tr in parents)
                    {
                        if (double.IsNaN(tr._flow))
                            tr._flow = 1.0;
                        jointInFlowSum += tr._flow;                            
                    }
                    foreach (TrackResult tr in children)
                    {
                        if (double.IsNaN(tr._flow))
                            unknownCount ++;
                    }
                    foreach (TrackResult tr in children)
                    {
                        if (unknownCount != 1)
                        {
                            if (double.IsNaN(tr._flow))
                            {
                                tr._flow = 1;
                                unknownCount--;
                            }
                            jointOutFlowSum += tr._flow;
                        }
                        else
                        {
                            tr._flow = jointInFlowSum - jointOutFlowSum;
                        }                        
                    }

                    // a joint
                    NetworkJoint joint = new NetworkJoint();
                    joint._children = new List<TrackResult>(children);
                    joint._parents = new List<TrackResult>(parents);

                    networkJointList.Add(joint);

                    downstreamSearch.Remove(current);

                    parents.Remove(current);
                    foreach (TrackResult tr in parents)
                    {
                        if (!tr._isRoot)
                        {
                            isIn = false;
                            foreach (TrackResult tr1 in upstreamSearch)
                            {
                                if (tr1._index == tr._index)
                                {
                                    isIn = true;
                                    break;
                                }
                            }
                            if (!isIn)
                                upstreamSearch.Add(tr);
                        }
                    }
                    foreach (TrackResult tr in children)
                    {
                        if (!tr._isLeaf)
                        {
                            isIn = false;
                            foreach (TrackResult tr1 in downstreamSearch)
                            {
                                if (tr1._index == tr._index)
                                {
                                    isIn = true;
                                    break;
                                }
                            }
                            if (!isIn)
                                downstreamSearch.Add(tr);   
                        }
                    }
                }

                // upstream search
                if (upstreamSearch.Count > 0)
                {
                    current = upstreamSearch[0];
                    currentIndex = current._index;
                    // get current's parents
                    parents.Clear();                    
                    foreach (TrackResult tr in trackResultList)
                    {
                        foreach (int rootIndex in current._rootIndexList)
                        {
                            if (rootIndex == tr._index && !parents.Contains(tr))
                            {
                                parents.Add(tr);
                            }
                        }
                    }

                    // get parents' children, i.e., current's siblings
                    children.Clear();
                    children.Add(current);                    
                    foreach (TrackResult parent in parents)
                    {
                        foreach (TrackResult tr in trackResultList)
                        {
                            foreach (int rootIndex in tr._rootIndexList)
                            {
                                if (rootIndex == tr._index && !children.Contains(tr))
                                {
                                    children.Add(tr);
                                    break;
                                }
                            }
                        }
                    }

                    // assign flow, upstream, => start from children
                    jointInFlowSum = 0.0;
                    jointOutFlowSum = 0.0;
                    unknownCount = 0;
                    foreach (TrackResult tr in children)
                    {
                        if (double.IsNaN(tr._flow))
                            tr._flow = 1.0;
                        jointOutFlowSum += tr._flow;
                    }
                    foreach (TrackResult tr in parents)
                    {
                        if (double.IsNaN(tr._flow))
                            unknownCount++;
                    }
                    foreach (TrackResult tr in parents)
                    {
                        if (unknownCount != 1)
                        {
                            if (double.IsNaN(tr._flow))
                            {
                                tr._flow = 1;
                                unknownCount--;
                            }
                            jointInFlowSum += tr._flow;
                        }
                        else
                        {
                            tr._flow = jointOutFlowSum - jointInFlowSum;
                        }
                    }

                    // a joint
                    NetworkJoint joint = new NetworkJoint();
                    joint._children = new List<TrackResult>(children);
                    joint._parents = new List<TrackResult>(parents);

                    networkJointList.Add(joint);

                    upstreamSearch.Remove(current);
                    children.Remove(current);
                    foreach (TrackResult tr in children)
                    {
                        if (!tr._isLeaf)
                        {
                            isIn = false;
                            foreach (TrackResult tr1 in downstreamSearch)
                            {
                                if (tr1._index == tr._index)
                                {
                                    isIn = true;
                                    break;
                                }
                            }

                            if (!isIn)
                                downstreamSearch.Add(tr);
                        }
                    }
                    foreach (TrackResult tr in parents)
                    {
                        if (!tr._isRoot)
                        {
                            isIn = false;
                            foreach (TrackResult tr1 in upstreamSearch)
                            {
                                if (tr1._index == tr._index)
                                {
                                    isIn = true;
                                    break;
                                }
                            }
                            if (!isIn)
                                upstreamSearch.Add(tr);
                        }
                    }
                }

                if (downstreamSearch.Count == 0 && upstreamSearch.Count == 0)
                    break;
            }

            return networkJointList;
        } // InitializeFlow()

        public void AssignK()
        {
            AssignK(_dstTrackResultList);
        } // AsignK()

        /// <summary>
        /// assign k = L / R ^ 4
        /// </summary>
        /// <param name="trackResultList"></param>
        private void AssignK(List<TrackResult> trackResultList)
        {
            foreach (TrackResult tr in trackResultList)
                tr._k = (double)tr._cpArray.Count / System.Math.Pow(tr._mean_r, 4);
        } // AssignK()

        public void FindRoots()
        {
            FindRoots(_dstTrackResultList, _srcImage);
        } // FindRoots()

        /// <summary>
        /// find root nodes
        /// </summary>
        /// <param name="trackResultList"></param>
        private void FindRoots(List<TrackResult> trackResultList, Bitmap srcImage)
        {
            int width = srcImage.Width;
            int height = srcImage.Height;
            int halfHeight = height / 2;
            double cpX, cpY;

            TrackResult current = new TrackResult();
            TrackResult parent = new TrackResult();
            TrackResult grandparent = new TrackResult();
            List<TrackResult> siblingList = new List<TrackResult>();

            TrackResult root;
            foreach (TrackResult tr in trackResultList)
            {
                if (tr._index == 1)
                    root = tr;
            }

            foreach (TrackResult tr in trackResultList)
            {
                siblingList.Clear();
                current = tr;
                if (current._rootIndexList[0] == -1)
                {
                    current._isRoot = true;
                    continue;
                }

                //if (current._rootIndexList[0] == 1)
                //{
                //    continue;
                //}

                if (!current._isLeaf)
                    continue;

                cpX = current._cpArray[0].x;
                cpY = current._cpArray[0].y;
                if (cpY < halfHeight)
                {
                    if (cpX > cpY || width - cpX > cpY)
                        current._isRoot = true;
                }

                cpX = current._cpArray[current._cpArray.Count - 1].x;
                cpY = current._cpArray[current._cpArray.Count - 1].y;
                if (cpY < halfHeight)
                {
                    if (cpX > cpY || width - cpX > cpY)
                        current._isRoot = true;
                }

                if (current._isRoot)
                {
                    // check current's siblings and get parent
                    foreach (TrackResult tr1 in trackResultList)
                    {
                        if (tr1._rootIndexList[0] == current._rootIndexList[0]
                            && tr1._index != current._index)
                            siblingList.Add(tr1);
                        if (tr1._index == current._rootIndexList[0])
                            parent = tr1;
                    }

                    if (siblingList.Count > 0)
                    {
                        foreach (TrackResult sibling in siblingList)
                            sibling._rootIndexList.Add(current._index);
                        current._rootIndexList.Clear();
                        current._rootIndexList.Add(-1);
                    }
                    else
                    {
                        // get grandparent
                        foreach (TrackResult tr2 in trackResultList)
                        {
                            if (tr2._index == parent._rootIndexList[0])
                            {
                                grandparent = tr2;
                                break;
                            }
                        }
                        // update parent's root index list
                        parent._rootIndexList.Clear();
                        for (int i = 1; i < current._rootIndexList.Count; i++)
                            parent._rootIndexList.Add(current._rootIndexList[i]);
                        parent._rootIndexList.Add(current._index);
                        // update current's root index list
                        current._rootIndexList.Clear();
                        current._rootIndexList.Add(-1);

                        // trace back
                        while (true)
                        {
                            current = parent;
                            parent = grandparent;
                            grandparent = null;

                            // check current's siblings, get grandparent meantime;
                            siblingList.Clear();
                            foreach (TrackResult tr3 in trackResultList)
                            {
                                if (tr3._rootIndexList[0] == parent._index
                                    && tr3._index != current._index)
                                    siblingList.Add(tr3);
                                if (tr3._index == parent._rootIndexList[0])
                                    grandparent = tr3;
                            }

                            if (siblingList.Count == 0)
                            {
                                // update parent root index list
                                parent._rootIndexList.Clear();
                                parent._rootIndexList.Add(current._index);
                            }
                            else
                            {
                                // update siblings root index list
                                foreach (TrackResult sibling in siblingList)
                                    sibling._rootIndexList.Add(current._index);
                                break;
                            }

                            if (grandparent == null)
                                break;
                        }
                    }
                } // if (current._isRoot)
            } // foreach (TrackResult current in trackResultList)
        } // FindRoots()

        public void FindLeaves()
        {
            FindLeaves(_dstTrackResultList);
        } // FindLeaves()

        /// <summary>
        /// find leave node
        /// </summary>
        /// <param name="trackResultList"></param>
        private void FindLeaves(List<TrackResult> trackResultList)
        {
            bool isLeaf;
            foreach (TrackResult tr in trackResultList)
            {
                isLeaf = true;
                foreach (TrackResult trChild in trackResultList)
                {
                    if (trChild._rootIndexList[0] == tr._index)
                    {
                        isLeaf = false;
                        break;
                    }
                }
                tr._isLeaf = isLeaf;
            }
        } // FindLeaves()

        public void ExtendBranchesToRoot()
        {
            ExtendBranchesToRoot(_dstTrackResultList, _srcImage);
        } // ExtendBranchesToRoot()

        /// <summary>
        /// extend branch starts to its root
        /// </summary>
        /// <param name="trackResultList"></param>
        /// <param name="srcImage"></param>
        private void ExtendBranchesToRoot(List<TrackResult> trackResultList, Bitmap srcImage)
        {
            // get threshold
            ConvertBitmapMatrix converter = new ConvertBitmapMatrix();
            Bitmap tempImage = new Bitmap(srcImage);
            eee.Sheffield.PZ.Imaging.Image.StandardizeImageTo8bbpIndexed(ref tempImage);
            PZMath_matrix srcMatrix = converter.ConvertGraylevelBitmapToMatrix(tempImage);
            tempImage.Dispose();
            SISThreshold sisThreshold = new SISThreshold();
            sisThreshold.ApplyDoubleGraylevelImageFilter(srcMatrix);
            int threshold = sisThreshold.Threshold;

            int rootIndex;
            TrackResult rootTR;
            PZPoint rootEndCP = new PZPoint();
            PZPoint branchStartCP;
            PZPoint candidateCP;
            PZDirection extendDirection;
            PZDirection candidateDirection = new PZDirection();
            int dX;
            int dY;
            int candidateX;
            int candidateY;
            int candidateIntensity;
            int currentIntensity;
            int minIntensity = int.MaxValue;
            bool toExtend;
            double distance;
            double angle;

            int width = srcImage.Width;
            int height = srcImage.Height;

            foreach (TrackResult tr in trackResultList)
            {
                rootIndex = tr._rootIndexList[0];
                if (rootIndex != -1)
                {
                    // search the root
                    foreach (TrackResult tr1 in trackResultList)
                    {
                        if (tr1._index == rootIndex)
                        {
                            rootTR = tr1;
                            // find the root end CP
                            rootEndCP = rootTR._cpArray[rootTR._cpArray.Count - 1];
                            break;
                        }
                    }
                   
                    while (true)
                    {
                        minIntensity = int.MaxValue;
                        toExtend = false;

                        // find the branch start CP
                        branchStartCP = tr._cpArray[0];
                        distance = branchStartCP.Distance(rootEndCP);
                        if (distance < 1)
                            break;

                        currentIntensity = (int)srcMatrix[(int)System.Math.Round(branchStartCP.y), (int)System.Math.Round(branchStartCP.x)];
                        
                        // direction extending branch to its root
                        extendDirection = new PZDirection(branchStartCP, rootEndCP);

                        // trial extension to 8 directions
                        #region try 8 directions
		                // (1, 1)
                        dX = 1;
                        dY = 1;
                        candidateX = (int)System.Math.Round(branchStartCP.x + dX);
                        candidateY = (int)System.Math.Round(branchStartCP.y + dY);
                        if (candidateX >= 0 && candidateX < width && candidateY >= 0 && candidateY < height)
                        {
                            candidateIntensity = (int)srcMatrix[candidateY, candidateX];
                            if (candidateIntensity < threshold)
                            {
                                angle = extendDirection.GetAngle(new PZDirection(dX, dY)) ;
                                if (angle < System.Math.PI / 2.0
                                    && candidateIntensity < minIntensity)
                                {
                                    minIntensity = candidateIntensity;
                                    candidateDirection = new PZDirection(dX, dY);
                                    toExtend = true;
                                }
                            }
                        }
                        // (1, 0)
                        dX = 1;
                        dY = 0;
                        candidateX = (int)System.Math.Round(branchStartCP.x + dX);
                        candidateY = (int)System.Math.Round(branchStartCP.y + dY);
                        if (candidateX >= 0 && candidateX < width && candidateY >= 0 && candidateY < height)
                        {
                            candidateIntensity = (int)srcMatrix[candidateY, candidateX];
                            if (candidateIntensity < threshold)
                            {
                                angle = extendDirection.GetAngle(new PZDirection(dX, dY));
                                if (angle < System.Math.PI / 2.0
                                    && candidateIntensity < minIntensity)
                                {
                                    minIntensity = candidateIntensity;
                                    candidateDirection = new PZDirection(dX, dY);
                                    toExtend = true;
                                }
                            }
                        }
                        // (1, -1)
                        dX = 1;
                        dY = -1;
                        candidateX = (int)System.Math.Round(branchStartCP.x + dX);
                        candidateY = (int)System.Math.Round(branchStartCP.y + dY);
                        if (candidateX >= 0 && candidateX < width && candidateY >= 0 && candidateY < height)
                        {
                            candidateIntensity = (int)srcMatrix[candidateY, candidateX];
                            if (candidateIntensity < threshold)
                            {
                                angle = extendDirection.GetAngle(new PZDirection(dX, dY));
                                if (angle < System.Math.PI / 2.0
                                    && candidateIntensity < minIntensity)
                                {
                                    minIntensity = candidateIntensity;
                                    candidateDirection = new PZDirection(dX, dY);
                                    toExtend = true;
                                }
                            }
                        }
                        // (0, 1)
                        dX = 0;
                        dY = 1;
                        candidateX = (int)System.Math.Round(branchStartCP.x + dX);
                        candidateY = (int)System.Math.Round(branchStartCP.y + dY);
                        if (candidateX >= 0 && candidateX < width && candidateY >= 0 && candidateY < height)
                        {
                            candidateIntensity = (int)srcMatrix[candidateY, candidateX];
                            if (candidateIntensity < threshold)
                            {
                                angle = extendDirection.GetAngle(new PZDirection(dX, dY));
                                if (angle < System.Math.PI / 2.0
                                    && candidateIntensity < minIntensity)
                                {
                                    minIntensity = candidateIntensity;
                                    candidateDirection = new PZDirection(dX, dY);
                                    toExtend = true;
                                }
                            }
                        }
                        // (0, -1)
                        dX = 1;
                        dY = -1;
                        candidateX = (int)System.Math.Round(branchStartCP.x + dX);
                        candidateY = (int)System.Math.Round(branchStartCP.y + dY);
                        if (candidateX >= 0 && candidateX < width && candidateY >= 0 && candidateY < height)
                        {
                            candidateIntensity = (int)srcMatrix[candidateY, candidateX];
                            if (candidateIntensity < threshold)
                            {
                                angle = extendDirection.GetAngle(new PZDirection(dX, dY));
                                if (angle < System.Math.PI / 2.0
                                    && candidateIntensity < minIntensity)
                                {
                                    minIntensity = candidateIntensity;
                                    candidateDirection = new PZDirection(dX, dY);
                                    toExtend = true;
                                }
                            }
                        }
                        // (-1, 1)
                        dX = -1;
                        dY = 1;
                        candidateX = (int)System.Math.Round(branchStartCP.x + dX);
                        candidateY = (int)System.Math.Round(branchStartCP.y + dY);
                        if (candidateX >= 0 && candidateX < width && candidateY >= 0 && candidateY < height)
                        {
                            candidateIntensity = (int)srcMatrix[candidateY, candidateX];
                            if (candidateIntensity < threshold)
                            {
                                angle = extendDirection.GetAngle(new PZDirection(dX, dY));
                                if (angle < System.Math.PI / 2.0
                                    && candidateIntensity < minIntensity)
                                {
                                    minIntensity = candidateIntensity;
                                    candidateDirection = new PZDirection(dX, dY);
                                    toExtend = true;
                                }
                            }
                        }
                        // (-1, 0)
                        dX = -1;
                        dY = 0;
                        candidateX = (int)System.Math.Round(branchStartCP.x + dX);
                        candidateY = (int)System.Math.Round(branchStartCP.y + dY);
                        if (candidateX >= 0 && candidateX < width && candidateY >= 0 && candidateY < height)
                        {
                            candidateIntensity = (int)srcMatrix[candidateY, candidateX];
                            if (candidateIntensity < threshold)
                            {
                                angle = extendDirection.GetAngle(new PZDirection(dX, dY));
                                if (angle < System.Math.PI / 2.0
                                    && candidateIntensity < minIntensity)
                                {
                                    minIntensity = candidateIntensity;
                                    candidateDirection = new PZDirection(dX, dY);
                                    toExtend = true;
                                }
                            }
                        }
                        // (-1, -1)
                        dX = -1;
                        dY = -1;
                        candidateX = (int)System.Math.Round(branchStartCP.x + dX);
                        candidateY = (int)System.Math.Round(branchStartCP.y + dY);
                        if (candidateX >= 0 && candidateX < width && candidateY >= 0 && candidateY < height)
                        {
                            candidateIntensity = (int)srcMatrix[candidateY, candidateX];
                            if (candidateIntensity < threshold)
                            {
                                angle = extendDirection.GetAngle(new PZDirection(dX, dY));
                                if (angle < System.Math.PI / 2.0
                                    && candidateIntensity < minIntensity)
                                {
                                    minIntensity = candidateIntensity;
                                    candidateDirection = new PZDirection(dX, dY);
                                    toExtend = true;
                                }
                            }
                        }
	#endregion

                        if (toExtend)
                        {
                            bool toAddCandidateCP = true;
                            candidateCP = new PZPoint(branchStartCP.x + candidateDirection.x, branchStartCP.y + candidateDirection.y);
                            // extended CP should be not in the track result
                            foreach (TrackResult tr2 in trackResultList)
                            {
                                if (tr2._index == tr._index)
                                    continue;
                                if (tr2.Contains(candidateCP))
                                {
                                    toAddCandidateCP = false;
                                    break;
                                }
                            }

                            // extended CP should shorten the distance to the root end CP
                            if (candidateCP.Distance(rootEndCP) > distance)
                                toAddCandidateCP = false;

                            if (toAddCandidateCP)
                            {
                                tr._cpArray.Insert(0, candidateCP);
                                tr._rArray.Insert(0, tr._mean_r);
                            }
                            else
                                break;
                        }
                        else
                            break;
                    } // while

                } // if (rootIndex != -1)
            } // foreach

        } // ExtendBranchesToRoot()

        /// <summary>
        /// remove duplicates
        /// </summary>
        public void RemoveDuplicates()
        {
            _dstTrackResultList = RemoveDuplicates(_srcTrackResultList);
        } // RemoveDuplicates()

        /// <summary>
        /// remove duplicates
        /// </summary>
        /// <param name="srcTrackResultList"></param>
        /// <returns></returns>
        private List<TrackResult> RemoveDuplicates(List<TrackResult> srcTrackResultList)
        {
            List<TrackResult> dstTrackResultList = new List<TrackResult>(srcTrackResultList);
            TrackResult containerTR = new TrackResult();
            bool isIn;
            for (int i = 0; i < dstTrackResultList.Count; i++)
            {
                isIn = false;
                TrackResult dstTR = dstTrackResultList[i];
                foreach (TrackResult srcTR in srcTrackResultList)
                {
                    if (!dstTR.Equals(srcTR))
                    {
                        if (dstTR.IsIn(srcTR))
                        {
                            containerTR = srcTR;
                            isIn = true;
                            break;
                        }
                    }
                }

                if (isIn)
                {
                    // remove dstTr from dstTrackResultList
                    dstTrackResultList.Remove(dstTR);

                    // update root index
                    foreach (TrackResult trTemp in srcTrackResultList)
                    {
                        if (trTemp._rootIndexList[0] == dstTR._index)
                        {
                            trTemp._rootIndexList[0] = dstTR._rootIndexList[0];
                            trTemp._layer--;
                        }
                    }

                    i--;
                }
            }
            return dstTrackResultList;
        } // RemoveDuplicates()

        public void RemoveShortResults()
        {
            RemoveShortResults(ref _dstTrackResultList, _lengthThreshold);
        } // RemoveShortResults()

        /// <summary>
        /// remove short results
        /// </summary>
        /// <param name="srcTrackResultList"></param>
        private void RemoveShortResults(ref List<TrackResult> srcTrackResultList, int lengthThreshold)
        {
            for (int i = 0; i < srcTrackResultList.Count; i++)
            {
                TrackResult dstTR = srcTrackResultList[i];
                if (dstTR._cpArray.Count <= lengthThreshold)
                {
                    srcTrackResultList.Remove(dstTR);
                    
                    // update root index
                    foreach (TrackResult trTemp in srcTrackResultList)
                    {
                        if (trTemp._rootIndexList[0] == dstTR._index)
                        {
                            trTemp._rootIndexList[0] = dstTR._rootIndexList[0];
                            trTemp._layer--;
                        }
                    }
                    
                    i--;
                }
            }
        } // RemoveShortResults()

        public void RemoveDisconnectedCPs(int distance)
        {
            RemoveDisconnectedCPs(ref _dstTrackResultList, distance);
        } // RemoveDisconnectedCPs()

        /// <summary>
        /// remove disconnected CPs from track results, given distance tolerance
        /// </summary>
        /// <param name="srcTrackResultIst"></param>
        /// <param name="distance"></param>
        private void RemoveDisconnectedCPs(ref List<TrackResult> srcTrackResultIst, int distance)
        {
            foreach (TrackResult tr in srcTrackResultIst)
            {
                int index = tr._cpArray.Count - 1;
                int indexM1 = index - 1;
                if (tr._cpArray[index].Distance(tr._cpArray[indexM1]) > distance)
                    tr.RemoveAt(index);
            }
        } // RemoveDisconnectedCPs()


        public void CheckOverlap(string fileName)
        {
            CheckOverlap(_dstTrackResultList, fileName, _srcImage);
        } // CheckLayer()

        /// <summary>
        /// check layer
        /// </summary>
        /// <param name="srcTrackResultList"></param>
        /// <param name="fileName"></param>
        private void CheckOverlap(List<TrackResult> srcTrackResultList, string fileName, Bitmap srcImage)
        {
            FileStream fs = new FileStream(fileName, FileMode.Create, FileAccess.Write);
            StreamWriter w = new StreamWriter(fs);

            List<TrackResult> temp;
            List<int> layerList = new List<int>(0);
            List<int> indexList = new List<int>(0);
            List<int> lengthList = new List<int>(0);
            bool findMatch;
            foreach (TrackResult tr in srcTrackResultList)
            {
                findMatch = false;
                foreach (PZPoint cp in tr._cpArray)
                {
                    layerList.Clear();
                    lengthList.Clear();
                    indexList.Clear();
                    temp = FindCPInTrackResults(cp, srcTrackResultList, 1.5);
                    
                    foreach (TrackResult tempTR in temp)
                    {
                        if (!indexList.Contains(tempTR._index))
                        {
                            indexList.Add(tempTR._index);
                            lengthList.Add(tempTR._cpArray.Count);
                        }
                    }
                    if (indexList.Count > 1)
                    {
                        findMatch = true;
                        // output
                        w.Write(cp.ToString() + ": index: ");
                        foreach (int i in indexList)
                        {
                            w.Write(i + " ");
                        }
                        w.WriteLine();

                        // draw overlapping segments
                        //string path = Path.GetDirectoryName(fileName);                        
                        //foreach (int i in indexList)
                        //{
                        //    foreach (TrackResult trSearch in srcTrackResultList)
                        //    {
                        //        if (trSearch._index == i)
                        //        {
                        //            Bitmap dstImage = new Bitmap(srcImage);
                        //            string imageFileName = path + "\\overlapping" + i + ".bmp";
                        //            foreach (PZPoint p in trSearch._cpArray)
                        //            {
                        //                int x = (int)System.Math.Round(p.x);
                        //                int y = (int)System.Math.Round(p.y);
                        //                dstImage.SetPixel(x, y, Color.Red);
                        //            }
                        //            dstImage.Save(imageFileName, ImageFormat.Bmp);
                        //        }
                        //    }
                        //}                                               
                    }
                }
                if (findMatch)
                {
                    w.Write("duplicate length: " + indexList.Count + " segment length: ");
                    foreach (int i in lengthList)
                        w.Write(i + " ");
                    w.WriteLine();
                }
            }

            w.Flush();
            w.Close();
            fs.Close();
        } // CheckLayer()


        public List<TrackResult> FindCPInTrackResults(PZPoint cp)
        {
            return FindCPInTrackResults(cp, _dstTrackResultList, 2);
        } // FindCPInTrackResults()

        /// <summary>
        /// find the track results to which a centre point belongs
        /// </summary>
        /// <param name="cp"></param>
        /// <param name="srcTrackResultList"></param>
        /// <returns></returns>
        private List<TrackResult> FindCPInTrackResults(PZPoint cp, List<TrackResult> srcTrackResultList, double distance)
        {
            List<TrackResult> dstTrackResults = new List<TrackResult>(0);
            foreach (TrackResult tr in srcTrackResultList)
            {
                foreach (PZPoint p in tr._cpArray)
                {
                    if (p.Distance(cp) < distance)
                    {
                        dstTrackResults.Add(tr);
                        break;
                    }
                }
            }
            return dstTrackResults;
        } // FindCPInTrackResults()


        public void IndexAll(List<TrackResult> trackResultList)
        {
            int i = 0;
            foreach (TrackResult tr in trackResultList)
            {
                tr._index = i;
                i++;
            }
        } // IndexAll()

        public int TotalLength()
        {
            return TotalLength(_dstTrackResultList);
        }

        private int TotalLength(List<TrackResult> trackResults)
        {
            int totalLength = 0;
            foreach (TrackResult tr in trackResults)
            {
                totalLength += tr._cpArray.Count;
            }
            return totalLength;
        }
        #endregion

        #region I/O method
        /// <summary>
        /// read track results from files
        /// </summary>
        public void ReadFile()
        {
            if (_pathName != "")
            {
                _srcTrackResultList = ReadFile(_pathName);
                _dstTrackResultList = new List<TrackResult>(_srcTrackResultList);
            }
        } // ReadFile()

        /// <summary>
        /// read track results from files
        /// </summary>
        /// <param name="pathName"></param>
        /// <returns></returns>
        private List<TrackResult> ReadFile(string pathName)
        {
            List<TrackResult> trackResultList = new List<TrackResult>();
            string[] fileNameArray = Directory.GetFiles(pathName, "*TrackResult*.txt");
            if (fileNameArray.Length > 0)
            {
                foreach (string fileName in fileNameArray)
                {
                    TrackResult tr = new TrackResult();
                    tr.Read(fileName);
                    trackResultList.Add(tr);
                }
            }

            //IndexAll(trackResultList);

            return trackResultList;
        } // ReadFile()


        public Bitmap DrawTrackResultsCentreLines()
        {
            return DrawTrackResultsCentreLines(_dstTrackResultList, _srcImage);
        } // DrawTrackResults()

        /// <summary>
        /// draw track result to dstImage
        /// </summary>
        /// <param name="trackResults"></param>
        /// <param name="dstImage"></param>
        /// <returns></returns>
        public Bitmap DrawTrackResultsCentreLines(List<TrackResult> trackResults, Bitmap srcImage)
        {
            Bitmap dstImage = new Bitmap(srcImage);
            foreach (TrackResult tr in trackResults)
            {
                foreach (PZPoint cp in tr._cpArray)
                {
                    int x = (int)System.Math.Round(cp.x);
                    int y = (int)System.Math.Round(cp.y);
                    dstImage.SetPixel(x, y, Color.Red);
                }
            }
            return dstImage;
        } // DrawTrackResults()

        /// <summary>
        /// index, root index, layer, average radius, length 
        /// </summary>
        /// <param name="fileName"></param>
        public void WriteDiagram(string fileName)
        {
            FileStream fs = new FileStream(fileName, FileMode.Create, FileAccess.Write);
            StreamWriter w = new StreamWriter(fs);

            foreach (TrackResult tr in _dstTrackResultList)
            {
                w.WriteLine(tr.ToString());
            }

            w.Flush();
            w.Close();
            fs.Close();
        }
        #endregion

    }
}
